Three-dimensional aspects of fluid flows in channels. I. Meniscus and Thin Film 

regimes 
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We study the forced displacement of a fluid-fluid interface in a three-dimensional channel formed 
by two parallel solid plates. Using a Lattice-Boltzmann method, we study situations in which a 
slip velocity arises from diffusion effects near the contact line. The difference between the slip and 
channel velocities determines whether the interface advances as a meniscus or a thin film of fluid is 
left adhered to the plates. We find that this effect is controlled by the capillary and Peclet numbers. 
We estimate the crossover from a meniscus to a thin film and find good agreement with numerical 
results. The penetration regime is examined in the steady state. We find that the occupation fraction 
of the advancing finger relative to the channel thickness is controlled by the capillary number and 
the viscosity contrast between the fluids. For high viscosity contrast, Lattice-Boltzmann results 
agree with previous results. For zero viscosity contrast, we observe remarkably narrow fingers. The 
shape of the finger is found to be universal. 



I. INTRODUCTION 

Advancing fronts in fluid systems involve the motion 
of a fluid-fluid interface, a surface that lives in a three- 
dimensional world, and which is often constrained by a 
solid boundary. A typical example is that of an interface 
moving in a channel [J 0, Q- 

Examples of advancing fronts in channels are imbibi- 
tion, a process in which a wetting fluid invades the chan- 
nel due to an uncompensated capillary pressure, and the 
viscous fingering process [3, y, IJ] , where a low- viscosity (or 
high-density) fluid penetrates a high- viscosity (or low- 
density) one. 

The problem of viscous fingering in a channel has been 
widely studied in the framework of Hclc-Shaw theory A 
Hele-Shaw cell is the two-dimensional limiting case of a 
very thin channel, where the equations of motion are av- 
eraged over the channel thickness. This reduces the inter- 
face to a line, the leading interface, that lives in the plane 
of the cell. As a consequence, the approximation discards 
any effects arising from the full three-dimensional struc- 
ture of the interface. 

Nonetheless, penetration in the gap of a Hele-Shaw 
cell is a fundamental three-dimensional effect that has 
important repercussions in the viscous fingering prob- 
lem. As theoretical studies have pointed out[f|, a thin 
film of viscous fluid left adhered to the cell plates as 
the front advances modifies the capillary pressure at the 
leading interface, thus altering the front morphology. 
This has been confirmed in experiments of steady vis- 
cous fingers @, where the presence of a thin film led to 
fingers not predicted by the two-dimensional theory. In 
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the following paper, we will address the role of the thin 
film in viscous fingers. 

A thin wetting film is not the only consequence of 
a three-dimensional intcrfacial structure. In the con- 
text of liquid films spreading over dry substrates @, [1], 
where a two-dimensional approximation is typically ap- 
plied, three-dimensional effects are also important. For 
instance, the stability of a spreading front depends on 
the wetting properties of the fluid. Experimentally, it 
has been observed Q that a crossover from unstable to 
stable fronts occurs when the dynamic contact angle ex- 
ceeds 7r/2, a situation which renders the velocity field 
within the film three-dimensional. 

The problem of a moving interface in a three- 
dimensional channel must take into account a dynamic 
contact line, the intersection point between the fluid- 
fluid interface and the channel walls. In classic fluid 
mechanics, a moving contact line violates the usual no- 
stick boun dary condition, leading to a divergent viscous 
dissipation[10j. Hence, contact line dynamics must con- 
sider a regularizing mechanism of the viscous dissipa- 
tion singularity. A slip velocity in the vicinity of a 
driven contact line arises naturally in diffuse interface 
models of binary fluids (Tl|, regularizing the singularity. 
These models consist of the usual Navier-Stokes equa- 
tions coupled to a convection-diffusion equation of an 
order parameter [T3|. Diffuse interface effects enter the 
force balance equations in the shape of order parameter 
gradients that play the role of a Young force. As a re- 
sult, the contact line slips over the solid boundary (Til. [l3| . 
Away from the contact line, order parameter gradients 
vanish and the stick boundary condition is recovered. 
The size of the diffusion region, Id, is then a measure 
of how strong slip is for a given system and is clearly an 
important parameter. This size was estimated by Briant 
and Yeomans[3|, who characterize Id for the case of an 
interface subjected to shearing walls. They focused on 
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the dependence of lr>(L in their notation) on the model 
parameters, finding a scaling relation that was verified 
numerically. 

Important implications arising from a relatively large 
or small slip velocity compared to the leading interface 
velocity in forced fronts can be foreseen. Whenever both 
velocities are comparable, the interface should maintain 
a meniscus shape. Conversely, as the slip velocity be- 
comes small compared to the channel velocity the inter- 
face shape should develop as a finger, leaving a thin film 
of fluid adhered to the walls of the channel. 

In this paper we study the penetration process across 
the channel thickness. We study the motion of the full 
three-dimensional interface between two viscous fluids 
when it is subjected to a gravitational body force. We 
treat the case of a strictly flat leading interface, focusing 
only in the details that pertain to the channel thickness. 
We work with symmetric fluids as well as with fluids of 
different densities or viscosities. 

We focus on two principal matters. We first describe 
how the contact line and leading interface velocities are 
related, and propose the mechanisms that determine the 
velocity ratio. We find that the velocity ratio is controlled 
by the force balance at the interface and by diffusion 
effects localized at the contact line. 

Secondly, we study the thin film that forms inevitably 
in the case of small slip. In that case the front decouples 
from the contact line leading to the growth of a finger, 
even when the interface is linearly stable to the Raylcigh- 
Taylor instability. We find that the fraction of occupa- 
tion of the thin film relative to the channel thickness 
is a function of the capillary number and the viscosity 
contrast between the fluids. The high viscosity contrast 
case is validated by comparing our results to the numer- 
ical work of Halpern and Gaver]15j which is consistent 
with previous results of Taylor [16j. For fluids with zero 
viscosity contrast, it turns out that the finger width has 
much lower values than for the high viscosity contrast 
case at fixed capillary number. 

The morphology of our fingers is very much alike to the 
Saffman- Taylor finger shape, a prediction of the Hclc- 
Shaw theory. Nevertheless, it is important to stress 
that although the case of a flat leading interface is two- 
dimensional, the equations of motion are not equivalent 
to those of the Saffman- Taylor problem. Therefore, pen- 
etration in the channel thickness cannot be attributed 
to the Saffman- Taylor instability. Likewise, the selection 
rule of the steady state, i.e., the actual dependence of 
the thin film thickness with the front velocity, cannot 
be mapped to the theoretical predictions of the viscous 
fingering theory. 

We will address these matters by means of numerical 
simulations of the mesoscopic equations of the system. 
To do so, we take advantage of a powerful integration al- 
gorithm in fluid dynamics, the Lattice-Boltzmann scheme 
for binary fluids. 

The paper is organized as follows. In Sec. [II] we 
present the equations that govern the system in the meso- 



scopic regime. In Sec. IIIII we briefly present the Lattice- 
Boltzmann algorithm for binary fluid flows. Sec. IIV Al 
is dedicated to simulation results of the forced interface, 
from which two steady state regimes are found; a non- 
penetrating regime, in which the interface advances as a 
meniscus, and a penetrating one, in which a single fin- 
ger emerges and achieves steady state. In Sec. IIV Bl we 
present a scaling argument of the equations of motion 
that leads to an estimate of the ratio between the slip and 
front velocities. Such argument explains the crossover 
from one regime to the other. In Sec. IIV CI we extend 
our results to fluids of different densities or viscosities. 
Sec. IIVDI is devoted to the steady state finger. Finally, 
in Sec. [V] we present the conclusions of this work. 



II. GOVERNING EQUATIONS 

We consider a channel formed by two solid plates par- 
allel to the xy plane, each of length L and infinite width, 
located at positions z = and z = b. Initially, two flu- 
ids fill the channel and are separated by a flat interface 
perpendicular to the solid walls, as shown in Fig. [TJ The 
equilibrium contact angle is hence Be = 7r/2. Contact 
lines are located at z = and z = b, while the leading 
interface is located at z = 6/2. 

To circumvent the complications of the sharp interface 
formulation, we introduce a mesoscopic variable, <j)(r); 
and order parameter which is constant in the bulk of 
each fluid and varies smoothly across a diffuse interfacial 
region. Within this approximation, the equilibrium state 
of the system is described by a Helmholtz free energy fl2T| 

T{pM = J dr(F(0,p) + |(V0) 2 ). 

The first term in the integrand is a volume contribution, 
given by V(<f>,p) = A<j) 2 /2 + Bcjfi/A + p/Zhxp. The p 
dependent term corresponds to an ideal gas contribution, 
while the 4> dependent terms allow for the coexistence of 
two phases. The presence of an interface is accounted 
for by the last term in the integrand, which penalizes 
spatial variations of the order parameter by a factor k. 
Minimization of T leads to the chemical potential, 
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FIG. 1: Schematic representation of the system. The leading 
interface and contact line positions are indicated. 
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and total pressure tensor (l7|. 



III. LATTICE BOLTZMANN METHOD 
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where S is the diagonal matrix. The pressure tensor has 
an ideal contribution given by P = -5, and an order 
parameter contribution. In equilibrium, the order pa- 
rameter profile for the flat interface (sketched in Fig. [T|) 
is (f>*(x,z) — — (j} e q tanh(x/£), where 4> eq = (-A/B) 1 / 2 is 
the bulk equilibrium value of the order parameter and 
£ = (— k/2A) 1 ^ 2 is the length scale of the interfacial re- 
gion; this profile leads to the difference between equilib- 
rium values A<fr = 2<p eq and the energy per unit area of 
the interface, a = (— 8kA 3 /9B 2 ) 1 ^ 2 . Since the interface 
is diffuse, a choice for the nominal interface position has 
to be made. We choose the level surface = 0. 

The divergence of the pressure tensor yields the force 
per unit volume that acts on the fluid: —VP — 0V/z. 
The first term is the pressure gradient, while the sec- 
ond arises from order parameter inhomogeneities. Con- 
sequently, the Navicr-Stokes equations are[l2j], 

pfdtfi+v-Vv) = —VP — 0V> + nV 2 v + pg, (1) 

where v is the fluid velocity, r\ is the fluid viscosity and 
g is the acceleration of gravity. 

The dynamics of the order parameter are described by 
a convection-diffusion equation, 
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V</> = MV>, 



(2) 



where M is a mobility. For small deviations from the 
equilibrium configuration, an expansion of the chemical 
potential in powers of <p — 4>* yields a first order diffusion 
coefficient D = M(A+ B(j> 2 q ), so the relative importance 
of the advective and diffusive terms can be estimated 
through a Peclet number, Pe = \v ■ W(f>\/\DW 2 cf>\. 

The system can be represented as a sheet of fluid in 
the xz plane with periodic boundary conditions applied 
in the y direction. This is equivalent to a channel of 
infinite width in the y direction with a flat leading in- 
terface. Stick boundary conditions are imposed at the 
walls, v(x, z = 0) = v(x, z = b) = 0, while no flow 
boundary conditions are imposed for the order param- 
eter, 4>v(x, z = 0) = 4>v(x, z = b) — 0. At both ends 
of the channel the flow is homogeneous. This is ensured 
by setting d x pv(x = 0,z) = d x pv(x = L, z) = and 
d x (f>v(x = 0, z) = d x (j)v(x = L,z) = 0. 

Contact line dynamics arise from the diffuse nature 
of the interface, which allows for slip in the interfacial 
region by a diffusive mechanism. The size over which 
slip takes place, Ijj, is a function of the fluid properties 
and has been estimated by Briant and Yeomans[14l|. who 
have given a scaling relation, lo ~ (r/^Af/A^ 2 ) 1 / 4 . 



We solve numerically Eqs. (Q]) and © by means of 
the Lattice-Boltzmann algorithm presented in Ref. [l7l ]. 
The dynamics are introduced by discretized Boltzmann 
equations of two distribution functions, 
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and 



g i (r + c l ,t+ 1) -gi(r,t) = — -(g 
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(4) 



In these equations, fi and gi are distribution functions, 
where the index i counts over the model velocity set. 
Space is discretized as a cubic lattice where nodes are 
joined by velocity vectors, C{. Space and time units in 
Eqs. ^ and (01 arc set to unity. Likewise, the density 
of the fluids is set to one. We use the D3Q15 velocity 
set, which consists of fifteen velocity vectors: six of mag- 
nitude 1 that correspond to nearest neighbors, eight of 
magnitude that correspond to third-nearest neigh- 
bors and one of zero magnitude that accounts for rest 
particles. In the D3Q15 model the speed of sound is 
c 8 = l/v3- In Eqs. j3]) and j4|, distribution functions are 
first relaxed to equilibrium values, represented by f^ q and 

gl q , with relaxation timescales 17 and T g . The term F? 
is related to the external forcing. Following the collision 
stage, distribution functions are propagated to neighbor- 
ing sites. 

Hydrodynamic variables are defined through moments 
of the fi and gi. The local density and order parame- 
ter are given by V_V / 4 = p and V_V g, L = <fi. The fluid 
momentum and order parameter current, are defined as 
Si f& = Pv an d Si di^i = 4>v- Local conservation 
of mass and momentum is enforced through the con- 
ditions_X)i /r = p, EiSf = 4>, J2 l f! 9 Ci = pv and 
YJiffi^i = 4>v. In equilibrium, the pressure tensor and 

chemical potential arc defined as V_V f? q CiCi = pvv + Pt 
and V\ g1 q CjC, = MpS + 4>vv. 

The equilibrium distribution functions and the forcing 
term are written as expansions in powers of v^M^, i.e., 



I A[ + 3v ■ Ci + -vv : CiC t 



G '. CiCv 



9 3=g 
- pu„ [ A 9 V + 3v ■ ^ + -vv : - -v 2 + G : CiC, 



and 



F =4u v [l 



2t/ 



/ • Ci(l +V-Ci) - v- f 



Here, v stands for the three possible magnitudes of the 
Ci set. Coefficient values are ujq = 2/9, u>i — 1/9 and 
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lu^ = 1/72; A J Q = 9/2 - 7/2TrP, A[ = A 1 ^ = 1/pTrP 



- 21/2M/Z, = 

— S), where 1 is 



and G s = 9/(2p)P - 3STrP; A 9 = 9/2- 
A 9 ^ = 3M[i/p and G 9 = 9/(2p)M/j,(l 
the unit matrix. 

Eqs. ([T]) and can be recovered as a Chapman- 
Enskog expansion of Eqs. (O and ((4]) [18] . The Lattice- 
Boltzmann scheme maps to the hydrodynamic model 
through the relaxation timcscalcs, i.e., i] ~ (2t/ — l)/6 
and M = (r g — 1/2)M, and through the body force 

/ = 99- 

Solid boundaries in the Lattice-Boltzmann method are 
implemented by means of the well known bounce-back 
rules [H, [nj. In the lattice nodes that touch the solid, 
the propagation scheme is modified so the distribution 
functions are reflected to the fluid rather than absorbed 
by the solid. As a consequence, a stick condition for 
the velocity is recovered approximately halfway from the 
fluid node to the solid node. 
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IV. RESULTS 



We study the process of penetration across the channel 
thickness in the presence of a dynamic contact line. As 
we have explained above, fingering is expected whenever 
the slip velocity is small compared to the leading inter- 
face velocity. In our model, slip is controlled by diffusion 
in the vicinity of the contact line. To measure the im- 
portance of diffusivity we use a typical definition of the 
Peclet number, Pe = Ub/D, where U is the velocity of 
the leading interface. The other relevant control param- 
eter is the capillary number, which follows from the ratio 
between viscous and capillary forces, Ca = rjU/a. We 
focus on flows governed by viscous and capillary forces. 
To enforce this situation we neglect the convective term 
in Eq. ([T]). To assure that we work on the low Mach num- 
ber regime, the fluid velocity is restricted to U < 0.01. 
For the case of small slip, we expect a thin film regime 
typical of experiments. We characterize this regime in 
terms of the finger width, viscosity contrast and capil- 
lary number. We compare our results with other studies 
from the literature. 



A. Effect of diffusivity, surface tension and 
viscosity 

We first consider two fluids with equal viscosities and 
densities. The size of the interface is set to £ = 0.57, 
which has been previously verified to give sufficiently ac- 
curate results for the variation of </> and its spatial deriva- 
tives across the interface [l?}- Starting from a flat in- 
terface configuration, we perform a set of five runs at 
fixed forcing, viscosity and surface tension. For each run 
we choose a different diffusion coefficient, which we fix 
through the mobility. In terms of dimcnsionlcss numbers 
this corresponds to fix Ca and vary Pe. Parameter val- 



FIG. 2: Interface evolution in the channel thickness direction 
for varying diffusion strength. Time interval between inter- 
faces is St ~ 2.17 in b/U units. The thick profile in each 
figure corresponds to the latest time. Meniscus regime: (a) 
D = 0.073 and (b) D = 0.049. Finger regime : (c) D = 0.024, 
(d) D = 0.012 and (e) D = 0.009. 



ues are U = 5 x 10~ 3 , r) = 10 _1 and a = 4.6 x 10~ 3 , 
where U is the expected leading interface velocity, calcu- 
lated as U = b 2 pg/(8r]). Channel dimensions are b = 23 
and L = 500. 

In Fig. [5] we show a time sequence of the interface po- 
sition for each run. In our simulations, v y = so a flat 
leading interface is located at z = 6/2. Sequences (a) 
and (b) correspond to runs with the highest diffusion 
coefficients. In both cases a steady meniscus is clearly 
observed. It is also appreciable that the meniscus in se- 
quence (a), corresponding to the highest diffusivity, is 
less curved than the meniscus in sequence (b) . The next 
three sequences, (c), (d) and (e), show an abrupt change 
in the interface configuration. Instead of a meniscus, we 
observe a penetrating structure that emerges from the 
center of the channel leaving a thin film of fluid adhered 
to the solid plates. The finger width in runs (c)-(e) is 
approximately 17 lattice spacings. For the size of the 
interface used, the order parameter saturates to its equi- 
librium value at the solid surface. Nonetheless, to rule 
out any effects associated to the size of the interface, we 
have verified that the finger width (relative to the channel 
thickness) does not depend on b, as we will sec bellow. 

All runs achieve a steady state in which the velocity 
of the leading interface is constant. This velocity is the 
same for runs (a) and (b) and due to mass conservation 
is slightly larger (a few percent) for runs (c), (d) and (e). 
The capillary number is not affected much by this effect, 
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TABLE I: Parameter values for ij and a varying runs, U ~ 

5 x lrr 3 , D = 7.5 x icr 2 . 

a = 4.6 x 10~ 3 T) = IP" 1 



rj shape a shape 

0.1 meniscus 0.0044 meniscus 

0.2 finger 0.0037 meniscus 

0.4 finger 0.0032 meniscus 

0.6 finger 0.0027 finger 



and we will take it as constant. The relevant effect is 
associated to the variation of diffusivity. 

The velocity of the contact line increases with increas- 
ing diffusivity, as can be deduced from the contact line 
position in sequences (c), (d) and (e). Nevertheless, the 
velocity of the leading interface and the width of the pen- 
etrating finger are the same for all three runs. This is a 
direct confirmation of the fact that contact line dynam- 
ics are decoupled from leading interface dynamics in the 
presence of a thin film, as proposed by Park and Homsy 
in Ref.gi. 

It is clear from these runs that the crossover for pen- 
etration is set by the difference between the leading in- 
terface velocity, U, and the slip velocity at the contact 
line, v s . For a meniscus, v s = U, while penetration oc- 
curs whenever v s < U. As v s depends on the strength of 
diffusivity, we can draw as a conclusion that penetration 
can be achieved by increasing Pe. 

We now explore the effect of capillarity on the dynam- 
ics of the interface. To do so we force the interface at fixed 
velocity, diffusivity and surface tension(resp. viscosity) 
while we vary the viscosity (resp. surface tension). As a 
consequence, Pe is fixed while Ca is varied. 

Results are summarized in Table |TJ The first column 
shows parameters for runs in which the viscosity is var- 
ied. We observe that penetration occurs as rj increases. 
The second column in Table [J shows results for varying 
surface tension. We observe that penetration occurs as 
a is decreased. We can conclude that capillarity plays a 
similar role as diffusivity, as penetration occurs for low 
values of Ca. 



B. The Onset of Penetration 

Our results suggest that the crossover from the menis- 
cus regime to the thin film regime is controlled at least 
by two mechanisms. On the one hand, viscous stresses 
deform the interface. As a result surface tension tends 
to restore the interface shape to its equilibrium value. 
On the other hand, advection causes order parameter 
gradients. As a consequence, diffusivity generates a slip 
velocity at the contact line. In this section we will see 
that the balance between these mechanisms is controlled 
by Pe and Ca. 

Let us write the force balance per unit volume of fluid 
in the frame of reference of the interface. We introduce 



orthogonal curvilinear coordinates, s, the arclength along 
the curve <p = 0, and u, the normal distance to a point 
on this curve. In terms of these coordinates the normal 
component of Eq. ([T]) (in absence of inertial terms) is 

d u P = -<f>d u (i + ijV 2 v n + pg n , (5) 

where the subscript n stands for the normal component 
and the subscript u denotes differentiation with respect 
to u. 

The force per unit area acting on the interfacial region 
is obtained by integrating ([5]) across the interface: 

AP = -<t{k° - «f ) + { V V 2 v n + P g n ) f, (6) 

where the term <t(k® — k^) arises from the integration of 
the chemical potential term [12], with and being 
the dynamic and equilibrium curvatures, which are pos- 
itive for a bump protruding in the x direction. We have 
assumed that neither of the last two terms in the right 
hand side vary appreciably across the interface. Eq. ([5]) 
should be interpreted as the usual Gibbs-Thomson con- 
dition plus a dynamic term proportional to £, which van- 
ishes either in equilibrium or in the sharp interface limit. 

We will now examine Eq. (J5]) in the vicinity of the con- 
tact line. The mcsoscopic nature of the interface gives 
rise to a finite size region where diffusion is important. 
This results in a slip velocity, v s , for the contact line. We 
now reproduce the scaling argument presented in Ref . [bij 
to obtain the diffusion size, Id, and consequently v s . We 
will subsequently compare the slip velocity to the lead- 
ing interface velocity in terms of Ca and Pe, which are 
parameters that can be linked to experiments. 

The slip velocity and the size of the diffusion region fix 
the magnitude of viscous dissipation in Eq. ©, 

AP ~ -*(k? k?) + fe + A9n) e (7) 

Since in the contact line region the time variation of the 
order parameter is d<j)/dt ~ v s A(f>/£,, the order parameter 
variation obeys 

A* DAcfr 
v s - ~ - ir . (8) 

Using Eq. (jHJ) to eliminate Id from Eq. (O we get 

The last term in this expression is order £, while the 
rest of terms are of order £°. The term in the left hand 
side is the excess pressure drop caused by the curvature 
difference, which is small in our simulations. Neglecting 
both the pressure gradient and the body force we extract 
the following scaling law for the slip velocity: 

2 ajnl - kZ)D 
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FIG. 3: Peclet and capillary numbers for meniscus and thin 
film regimes. □ Meniscus regime. ■ Thin film regime, o 
Meniscus regime (different viscosities). • Thin film regime 
(different viscosities). A Meniscus regime (different densi- 
ties). A Thin film regime (different densities). The solid line 
corresponds to Eq. ([9}, with a ~ 0.3. 

The interface curvature is a consequence of the under- 
lying velocity profile, which is set by the thickness of 
the channel. Therefore, the curvature difference scales as 
(k®-K%) ~ a/b, with a being a typical amplitude. Using 
this expression and measuring v s in units of the leading 
interface velocity we arrive at the following expression: 

Q^) 2 „ aCa-iPe- 1 . 

This indicates that both Pe and Ca control how the slip 
velocity compares to the leading interface velocity. For a 
meniscus v s — U, so we arrive at the following condition: 

Pe = aCar x . (9) 

To check the validity of the prediction we perform sev- 
eral runs of forced interfaces varying simulation param- 
eters in a wide range(see Tables Hill and HV|) . We cover 
up to four decades in the Pe and Ca until numerical sta- 
bility issues of the code begin to show up. In Fig. [3] we 
show a plot of our results in the Ca~ 1 Pe plane. Data 
shown in the figure sketches two regions; at high Pe and 
Ca values the thin film regime is observed, whereas the 
meniscus corresponds to low values of both parameters. 
We also show our prediction, for which we find a fitting 
value for a, namely a ~ 0.3. 

Typical experiments with molecular liquids correspond 
to high, O(10 2 ) - O(10 3 ), values of the Peclet number [1, 
l2pl |. thus falling in the thin film region sketched in the 
diagram. However, menisci are expected in systems with 
a diffuse interface, such as colloid-polymer mixtures[2l|. 
In terms of experimental parameters, the diffusion length 
scales as Id ~ { T 1^ 2 ^/ cr ( K a ~~ K a )) 1 ^ 4 - In colloid-polymer 
mixtures the ratio (£ 2 /< 7 ) 1 ^ 4 is about 10 2 times larger 
than for molecular liquids. As a consequence, in such 
systems menisci should be observable at relatively high 
capillary numbers. For molecular liquids, this effect can 
be achieved by decreasing the system size, for instance, 



in microfluidic devices, where the typical size of the chan- 
nel is a few microns [23|. For such small sizes, the Peclet 
number is O(l), about 10 3 times smaller than for tradi- 
tional channels. Hence, the transition from menisci to 
thin films would be observable in the regime of relatively 
high capillary number, say Ca O(10 _1 ) — O(10°). 

C. Asymmetric Fluids 

We have shown that the ratio between the leading in- 
terface and contact line velocities is controlled by the in- 
terplay between the Peclet and capillary numbers. We 
expect that this fact holds for fluids of either differ- 
ent densities or viscosities. A forced interface between 
asymmetric fluids can be destabilized by virtue of the 
Rayleigh- Taylor instability, when the more dense fluid 
displaces the less dense one. We explore situations for 
which the instability is absent. To ensure this, we keep 
the channel thickness below the first unstable wave- 
length, given by l c = 2tt(lt/ Apg) 1 / 2 ^. 

The forcing is set according to / = 
8r]((j))U eX p/b 2 A((f))x. where the local viscos- 

ity is set according to the mixing rule r)(<t>) = 
(r)2 + ?7i)(l — C(j}/4> e q)/2, characterized by the vis- 
cosity contrast c = (i]2 — ?7i)/(?72 + ?7i), and U exp is 
the maximum expected velocity for a Poiseuille profile. 
The 4> dependent part is set as A((f>) = 1 if c ^ and 
A = ((f) + 4> eq ) / A(j) otherwise. In experiments the typical 
situation is that an effectively inviscid fluid displaces a 
viscous one, which corresponds to c — > 1. We approach 
this situation by setting c = 0.9. Following the general 
convention in the literature, here we define the capillary 
number as Ca = ^U/a. For all cases, we consider that 
both fluids have the same diffusion coefficient. 

We have performed a set of runs in which we vary 
Pe at fixed Ca for fluids with finite density or viscosity 
contrast. The details of the runs are summarized in Ta- 
bic |TTJ For each case, both menisci and fingers can be 
obtained depending on the value of the Peclet number. 
As expected, penetration occurs for sufficiently high Pe. 
We can conclude that the appearance of a thin film is 
independent of the Rayleigh- Taylor instability. In Fig. [3] 
we plot results of this section in the PeCa^ 1 plane. For 
fluids of different densities this value is consistent with 
the symmetric estimate of a ~ 0.3. For fluids of different 
viscosities the crossover occurs at a ~ 0.5. Anyhow, the 
qualitative behavior remains independent of the degree 
of asymmetry between the fluids. 

D. Steady state finger in the channel thickness 
direction 

We now turn our attention to the steady state finger 
that appears for high values of the product CaPe, which 
is the usual situation in most experiments. As we have 
shown in Sec. IIV Al diffusion only affects the motion of 
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TABLE II: Parameter values for runs of asymmetric fluids. 

Varying viscosities 
a = 4.6 x 10~ 3 , U ~ 2 x 10" 



Varying densities 

a = 9.2 x 10~ 3 , t/ ~ 4 



x 10 



-3 



D shape 
0.146 meniscus 
0.098 meniscus 
0.049 finger 
0.024 finger 
0.018 finger 



D 
0.0488 
0.0244 
0.0122 
0.006 
0.001 



shape 
meniscus 
meniscus 
meniscus 
meniscus 
finger 



the contact line, and has a negligible effect in the steady 
state finger. The relevant control parameters, as pro- 
posed in the literature, are the capillary number and the 
viscosity contrast between the fluids. The steady state 
if often characterized by measuring the finger width, Xb, 
which is the fraction of occupation of fluid 1 relative to 
the channel thickness. 

We explore A& as a function of Ca at a given c value. 
We consider three situations, c = and zero density 
contrast, c = with finite density contrast and c ^ 
with zero density contrast. For the last case we choose 
c = 0.90 and c = 0.95. The low Ca runs have been per- 
formed varying b to rule out lattice artifacts. We have 
found that results do not depend on the channel thick- 
ness chosen, the smallest thickness considered here being 
b = 23. Tabs. fVl IVII and lVIII displav the parameter values 
used in each run. 

In Fig. [4] we plot Ab as a function of Ca. We find 
that Xb depends on the viscosity contrast, as the c = 
points fall in a clearly different curve than the c = 0.9 
and c = 0.95 points. 

On the contrary, the density contrast does not play a 
relevant role. Points belonging to the c = curve were 
obtained using five different gap sizes; b = 147, b = 23, 
b = 35, b = 51 and b = 95, of which the b = 35 set 
was done with fluids of different densities. Results show 
no difference between zero or finite density difference. In 
fact, the gap size for the latter was large enough for the 



1 

0.9 
0.8 
Xb 0.7 
0.6 
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FIG. 4: Finger width as a function of the Ca. The error in the 
measured finger width is calculated from one lattice spacing 
and corresponds to approximately <5Af, ~ 0.04 in the figure. 



Ray lcigh- Taylor instability to be present. This means 
that the finger can develop as a consequence of low diffu- 
sion or by virtue of the Rayleigh- Taylor instability. Still, 
the steady state remains insensitive to the mechanism 
that leads to penetration and is selected by Ca and c. 

Previous analytic predictions correspond to the low 
Ca regime at c = 1 and were carried out first by 
Bretherton[24], who found that the finger width decays 
as At, — ► 1 — 1.337Ca 2 / 3 , as Ca — > 0. An extension 
was done by Taylor[l6j], up to Ca < 0.09, for which 
he reported a decaying exponent of one-half. Numeri- 
cally, Reinclt and Saffman (25| solved the Stokes equa- 
tions using a finite difference algorithm, and considered 
values up to Ca < 2, which match the one-half expo- 
nent of Taylor at small Ca. Halpern and Gaver[l5| ex- 
tended the prediction beyond Ca = 2 by means of a 
boundary element analysis of the Stokes equations. Their 
results can be fitted to an exponential law Xb = 1 — 
0.417 (l - cxp(-1.69Ca - 5025 )) (shown in Fig. |]§, which 
reproduces Rcinelt and Saffman results and matches the 
power law prediction of Taylor for low Ca. For large Ca, 
this law saturates to a limiting value of Ah = 0.583. Pre- 
vious Lattice-Boltzmann studies have also addressed this 
problem. Kang et a/[26[ studied the range 0.2 < Ca < 2, 
obtaining good agreement with Halpern and Gaver re- 
sults. Langaas and Yeomans[27j considered the range 
0.079 < Ca < 4.6 and were able to reproduce the results 
of Halpern and Gaver for Ca < 2. For Ca > 2 they ob- 
tained smaller finger widths than those of Halpern and 
Gaver. 

Our results cover up to five decades in the capillary 
number, 10~ 2 < Ca < 10 2 ^. They match Halpern 
and Gaver prediction as c — > 1. Already at c = 0.90, we 
reproduce accurately the finger width saturation value, 
for which we find A b = 0.573 ± 0.022. For small Ca, the 
error increases for the c = 0.9 runs. We improve this 
situation by increasing the viscosity contrast, as can be 
appreciated in Fig.[4j At Ca ~ 0.09 the error for c = 0.9 
is 4%,while for c = 0.95 it reduces to 2%. For Ca = 0.008 
the error is 2% at c = 0.95. This agreement shows that 
the Lattice-Boltzmann approach gives accurate results 
for a wide range of Ca, improving previous results [27|. 

As can be seen from Fig. 21 fingers with zero viscosity 
contrast, a case that has not been studied previously, 
are much narrow than fingers with c = 0.9 or c = 0.95. 
The dependence of the finger width with Ca has a power 
law behavior for 0.1 < Ca < 1, with an exponent m = 
0.29 ± 0.02. For Ca O(10), the finger width saturates to 
a notably small value, Ab ~ 0.386 ±0.014, which remains 
an open question. 

We now focus on the shape of the steady finger. Finger 
profiles shown in Fig. [5] are strongly reminiscent of the 
single finger solution of the Saffman- Taylor problem 
in which fingering occurs in the xy plane of a Helc-Shaw 
cell as a consequence of viscosity or density asymmetries 
between the fluids. Our problem is fundamentally dif- 
ferent because penetration in the channel thickness can 
occur for linearly stable interfaces in the context of hy- 
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drodynamic stability, e.g., by virtue of low diffusivity at 
the contact line. Moreover, even in the case where the in- 
terface is linearly unstable, it is due to a Rayleigh- Taylor 
instability and not through the Saffman- Taylor one. 

Still, we compare our finger profiles with the Saffman- 
Taylor ones. To do so, we recall the semi-empiric shape 
found by Pitts (29[, which for our geometry reads 



cos 



irz 



exp 



TTX 



(10) 



where x' and z' measure the distance from the finger tip. 
For both c = and c = 0.9, Eq. (flQ|) is a good approx- 
imation to our profiles at low Ca. This agreement is 
lost gradually as Ca increases. In Fig. [5] we show inter- 
face profiles for c = 0.0 and c = 0.9 at the smallest and 
largest Ca considered. Profiles for all other Ca values lie 
between the shown profiles and are omitted from the fig- 
ure. For c ~ 0.0, Eq. fTI)]) describes better our profiles for 
large Ca, while for c = 0.9 deviations from Pitts result 
arc observed as Ca increases. 



V. CONCLUSIONS 

We have studied the forced motion of a fluid-fluid inter- 
face in a three-dimensional channel by means of a meso- 
scopic model that takes into account contact line dynam- 
ics. 

Our results describe two possible scenarios regarding 
interface dynamics. A meniscus regime is found when- 
ever the contact line velocity is comparable to the lead- 
ing interface velocity. Conversely, when the contact line 
velocity is smaller than the leading interface velocity the 
meniscus configuration is lost, leading to penetration of 
one fluid into the other in a fingering fashion. A thin 




FIG. 5: Rescaled interface profiles for the lowest and highest 
Ca values of the c = 0.0 and c = 0.9 runs. For c = 0.0: 
oCa = 0.01 and »Ca = 74.68. For c = 0.9: UCa = 0.17 and 
■ : Ca = 52.82. Solid Line: Pitts semi-empirical finger shape. 
The error bar is shown in the bottom-right and corresponds 
to one lattice spacing. The length of the diffuse interface 
corresponds to approximately one half of a unit in the figure. 



film of displaced fluid is hence left adhered to the plates 
of the channel. 

The crossover from the meniscus to the thin film regime 
is controlled by the competition between surface and vis- 
cous stresses, as well as by the competition of diffusive 
and advective timescales, on top of the usual hydrody- 
namic instabilities. These mechanisms can be accounted 
for through simple scaling arguments. We find a pre- 
diction for the crossover in terms of the capillary and 
Peclet numbers which describes accurately our numeri- 
cal results. Menisci are found for low capillary and Peclet 
numbers, when surface tension and diffusion dominate 
over viscous stresses and advection respectively. An ex- 
ample of a system where diffusion is important is that 
of colloid-polymer mixtures. For such systems, the rel- 
atively large size of the interface together with low sur- 
face tensions leads to large diffusion regions near de con- 
tact line. For instance, in Ref. [2l[ the size of the in- 
terface is typically £ ~ 10 fim while the surface tension 
is (J — lfiN. For molecular liquids the size of the inter- 
face is of the order of nanometers, while a ~ 10 mN/m. 
As explained in Sec lIV Bl the size of the diffusion re- 
gion scales as Id ~ (£ 2 /ct) 1 / 4 , all other parameters kept 
constant. With these values the ratio of the diffusion 
length between colloid-polymer mixtures and molecular 
liquids is at least two orders of magnitude. Thin films 
are obtained for high values of both parameters, when 
advection and viscous stresses are dominant. Our pre- 
diction works for both symmetric and asymmetric flu- 
ids. From the crossover prediction, we propose that thin 
films can be assured as long as CaPe > 0.3 for symmet- 
ric fluids and CaPe > 0.5 for asymmetric fluids. Any- 
how, this values are consistent with the typical experi- 
mental regime. For instance, experiments of Ref. Q were 
done in a channel of thickness b = 7.95 x 10 -4 m, with 
a silicone oil with rj = 9.3 cP, a = 20.1 dyn/cm and 
D = 1.4946 x 10~ 6 cm 2 /s[13. Typical velocities in the 
experiments ranged from U = 0.01 m/s to U — 0.1 m/s. 
With these values, we obtain CaPe ~ 2.6 x 10 4 , which is 
consistent with our prediction for the thin film regime. 

We have examined the steady state of the thin film 
regime. We have found, in agreement with Ref. that 
contact line dynamics does not affect the steady state 
finger shape. The capillary number and the viscosity 
contrast between the fluids determine the shape of the 
finger. 

For a low- viscosity fluid pushing a high- viscosity one, 
the finger narrows with increasing capillary number down 
to a limiting value of 0.57 in units of the channel thick- 
ness. These results agree with the numerical results of 
Ref. [lH for the wide range of capillary numbers consid- 
ered. Due to computational limitations we do not inves- 
tigate fingers with very low capillary numbers. Nonethe- 
less, as we recover results from Ref. 15|, we expect that 
the low capillary number limit can be recovered by our 
method as well. 

For fluids with equal viscosities we have found a curve 
of the finger width as a function of the capillary number 
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that does not follow any previous results. The width of 
the finger decreases with increasing capillary number, an 
expected observation, but to a remarkably limiting width 
of 0.38 in units of the channel thickness. This contrasts 
with the saturation value of the asymmetric case. 

The steady state is independent of whether or not the 
interface is linearly unstable to the Raylcigh- Taylor insta- 
bility. This reinforces the conjecture of the steady state 
being independent of the mechanism that first leads to 
penetration of one fluid to the other. For low capillary 
numbers the shape of our fingers is universal and is con- 
sistent with the finger shape of Pitts, which suggests that 
the steady state can be described on simple mechanical 
equilibrium grounds. 

In a future work we will address the problem of vis- 
cous fingering allowing for a non-flat leading interface. 
As we have shown, the shape of the interface across the 
channel thickness can be controlled by tuning the diffu- 
sion strength in the contact line. Hence, it is possible to 
describe situations in which the leading interface under- 
goes a fingering process both in the presence and absence 
of a thin film in the channel thickness direction. In the 
presence of a thin film, additional control over the in- 
terface shape can be gained by choosing between fluids 
of equal or different viscosities. These features are very 
convenient to study in detail the three-dimensional effects 
that arise in the viscous fingering problem. 



TABLE III: Parameter values which develop a meniscus. For 
all runs b = 23. 



a 


D 


V 


u 


Ca 


Pe 


0.0092 


0.0976 


0.100 


0.0067 


0.0723 


1.56 


0.0092 


0.0488 


0.100 


0.0067 


0.0724 


3.12 


0.0092 


0.0976 


0.100 


0.0022 


0.0239 


0.52 


0.0092 


0.0488 


0.100 


0.0022 


0.0239 


1.04 


0.0092 


0.0244 


0.100 


0.0021 


0.0232 


2.02 


0.0092 


0.0122 


0.100 


0.0021 


0.0229 


3.98 


0.0092 


0.0244 


0.100 


0.0021 


0.0232 


2.00 


0.0092 


0.0122 


0.100 


0.0021 


0.0229 


3.98 


0.0092 


0.0066 


0.100 


0.0021 


0.0227 


7.86 


0.0092 


0.0031 


0.100 


0.0021 


0.0223 


15.5 


0.0092 


0.0976 


0.015 


0.0063 


0.0102 


1.48 


0.0092 


0.0488 


0.015 


0.0063 


0.0103 


2.97 


0.0092 


0.0244 


0.015 


0.0063 


0.0103 


5.94 


0.0092 


0.0122 


0.015 


0.0063 


0.0103 


11.9 


0.0092 


0.0091 


0.015 


0.0063 


0.0103 


15.8 


0.0092 


0.0091 


0.015 


0.0010 


0.0016 


2.52 


0.0009 


0.0240 


5.000 


0.0002 


1.0963 


0.20 


u.uuuy 




u.ouu 


n nnn9 
u.uuuz 


u.iuyo 


u.zu 


0.0009 


0.0240 


0.050 


0.0002 


0.0110 


0.20 


0.0044 


0.0750 


0.100 


0.0033 


0.0755 


1.02 


0.0037 


0.0750 


0.100 


0.0033 


0.0889 


1.02 


0.0032 


0.0750 


0.100 


0.0033 


0.1046 


1.02 


0.0027 


0.0750 


0.100 


0.0033 


0.1230 


1.02 


0.0044 


0.0750 


0.100 


0.0043 


0.0984 


1.32 


0.0037 


0.0750 


0.100 


0.0043 


0.1158 


1.32 
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TABLE IV: Parameter values which develop a finger. For all 
runs b = 23. 



a 


D 


V 


u 


Ca 


Pe 


0.0092 


0.0244 


0.1 


0.0080 


0.09 


7.54 


0.0092 


0.0122 


0.1 


0.0080 


0.09 


15.08 


0.0092 


0.0092 


0.1 


0.0080 


0.09 


20.10 


0.0092 


0.0092 


0.1 


0.0040 


0.04 


10.06 


0.0092 


0.0015 


0.1 


0.0030 


0.03 


45.24 


0.0092 


0.0976 


0.4 


0.0067 


0.29 


1.56 


0.0092 


0.0488 


0.4 


0.0067 


0.29 


3.16 


0.0092 


0.0244 


0.4 


0.0067 


0.29 


6.32 


0.0092 


0.0122 


0.4 


0.0067 


0.29 


12.64 


0.0092 


0.0092 


0.4 


0.0067 


0.29 


16.84 


0.0046 


0.0122 


0.1 


0.0050 


0.11 


4.72 


0.0046 


0.0061 


0.1 


0.0050 


0.11 


18.86 


0.0046 


0.0732 


0.2 


0.0050 


0.22 


1.58 


0.0046 


0.0732 


0.4 


0.0050 


0.43 


1.58 


0.0046 


0.0732 


0.6 


0.0050 


0.65 


1.58 


0.0046 


0.0732 


0.8 


0.0050 


0.87 


1.58 


0.0005 


0.0073 


0.8 


0.0050 


8.70 


15.72 


0.0000 


0.0008 


0.8 


0.0050 


86.98 


157.10 


0.0011 


0.0183 


0.8 


0.0125 


8.70 


15.72 


0.0000 


0.0008 


0.8 


0.0050 


86.98 


157.10 


0.0009 


0.0240 


1.0 


0.0010 


1.10 


0.96 


0.0009 


0.0192 


1.0 


0.0010 


1.10 


1.20 


0.0009 


0.0144 


1.0 


0.0010 


1.10 


1.60 


0.0009 


0.0120 


1.0 


0.0010 


1.10 


1.92 


0.0091 


0.0960 


1.0 


0.0100 


1.10 


2.40 


0.0009 


0.0060 


5.0 


0.0002 


1.10 


0.76 


0.0009 


0.0012 


5.0 


0.0002 


1.10 


3.84 


0.0027 


0.0750 


0.1 


0.0048 


0.18 


1.48 



TABLE VI: Parameter values for stationary fingers presented 



in Sec. ITVDl For all runs b = 35. 



a 


V 


U 


o a 


U 


1 e 




Non-Uniform Forcing, c 


= u.u 








b = 35 














0.0046 


0.50 


0.0049 


O.ooo 


U.000 1 1 


O A 1 1 *7 

241.1 ( 


U.OOZ 


0.0023 


0.50 


0.0051 


1.109 


0.00071 


250.09 


0.493 


0.0011 


0.50 


0.0052 


2.270 


0.00071 


255.98 


0.450 


0.0006 


0.50 


0.0053 


4.604 


0.00071 


259.57 


0.423 


0.0003 


0.50 


0.0053 


9.272 


0.00071 


261.38 


0.404 


0.00014 


0.50 


0.0054 


18.617 


0.00071 


262.39 


0.393 


0.00007 


0.50 


0.0054 


37.306 


0.00071 


262.91 


0.388 


0.00002 


0.50 


0.0016 


43.688 


0.00037 


154.85 


0.383 


0.00004 


0.50 


0.0054 


74.685 


0.00071 


263.16 


0.386 



TABLE V: Parameter values for stationary fingers presented 



in Sec. ITVDl 



a 


V 


V 


Ca 


D 


Pe 


A 6 


Uniform Forcing 


;, c = 0.0 










b = 23 














0.0044 


0.10 


0.0084 


0.192 


0.00185 


104.2 


0.661 


0.0037 


0.10 


0.0085 


0.228 


0.00222 


87.87 


0.650 


0.0032 


0.10 


0.0085 


0.270 


0.00265 


73.95 


0.641 


0.0027 


0.20 


0.0089 


0.662 


0.00064 


321.1 


0.561 


0.0023 


0.20 


0.0090 


0.790 


0.00076 


272.34 


0.509 


0.0019 


0.20 


0.0091 


0.936 


0.00091 


229.33 


0.490 


b = 51 














0.0044 


0.10 


0.0085 


0.194 


0.00442 


97.68 


0.678 


0.0037 


0.10 


0.0085 


0.230 


0.00449 


96.84 


0.653 


0.0032 


0.10 


0.0086 


0.272 


0.00459 


95.59 


0.640 


0.0027 


0.20 


0.0087 


0.646 


0.00490 


90.31 


0.532 


0.0023 


0.20 


0.0090 


0.789 


0.00489 


93.75 


0.524 


0.0019 


0.20 


0.0091 


0.941 


0.00490 


95.05 


0.497 


& = 95 














0.0044 


0.10 


0.0084 


0.192 


0.01608 


49.64 


0.706 


0.0037 


0.10 


0.0085 


0.229 


0.01621 


49.75 


0.682 


0.0032 


0.10 


0.0086 


0.272 


0.01634 


49.85 


0.661 


0.0027 


0.10 


0.0087 


0.323 


0.01648 


49.93 


0.658 


0.0023 


0.10 


0.0087 


0.383 


0.01661 


50 


0.622 


0.0019 


0.10 


0.0088 


0.455 


0.01675 


50.04 


0.598 


b = 147 














0.0046 


0.005 


0.0095 


0.010 


0.1464 


9.54 


0.835 


0.0046 


0.01 


0.0094 


0.020 


0.1464 


9.54 


0.828 



TABLE VII: Parameter values for stationary fingers presented 
in Sec. lIVDl 



a 


m 


U 


Ca 


D 


Pe 




Uniform Forcing 


c = 0.9 










b = 23 














0.0046 


0.38 


0.0092 


0.756 


0.00976 


21.57 


0.670 


0.0011 


0.38 


0.0097 


3.218 


0.00244 


91.76 


0.603 


0.0023 


0.38 


0.0095 


1.569 


0.00488 


44.73 


0.644 


0.0006 


0.38 


0.0099 


6.519 


0.00122 


185.91 


0.585 


0.0003 


0.38 


0.0099 


13.124 


0.00061 


374.28 


0.578 


0.0001 


0.38 


0.0100 


26.373 


0.00031 


752.09 


0.574 


0.0001 


0.38 


0.0100 


52.817 


0.00015 


1506.22 


0.572 


b = 51 


0.0046 


0.095 


0.0040 


0.082 


0.0061 


33.3 


0.805 


0.0046 


0.095 


0.0108 


0.223 


0.0061 


90.3 


0.743 


Uniform Forcing 


c = 0.95 








b = 147 














0.0046 


0.095 


0.0041 


0.084 


0.14 


4.07 


0.822 


0.0046 


0.049 


0.0007 


0.008 


0.14 


0.71 


0.931 


0.0046 


0.488 


0.0008 


0.087 


0.01 


9.89 


0.810 
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